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WAVE-PLAN ANALYSIS OF UNSTEADY FLOW 

by Don J. Wood, Robert G. Dorsch, 

and Charlene Lightner 

Lewis Research Center 
Cleveland, Ohio 

ABSTRACT: A digital distributed parameter model for computing unsteady 
flow in liquid-filled fluid systems is presented. The wave-plan method em- 
ployed in the model involves essentially the synthesis of incremental pressure 
pulses. The analysis is presented in a form general enough to be applied to a 
variety of hydraulic systems. To illustrate the application of the method to a 
specific system, the response of a long, straight liquid-filled line to sinus- 
oidal inlet flow and pressure perturbations is computed. Both a constant- 
cross-section line and a tapered line are analyzed in the example. 

KEY WORDS: hydraulics, unsteady flow, distributed parameter, nonlinear, 
digital. 
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WAVE-PLAN ANALYSIS OF UNSTEADY FLOW 


1 2 
by Don J. Wood , A. M. , A.S.C.E., Robert G. Dorsch , 

and Charlene Lightner^ 

SYNOPSIS 

An analytical method for computing unsteady flow conditions in liquid- 
filled flow systems is developed. The method which is called the wave plan in- 
corporates distributed parameter and nonlinear effects including the effects of 
viscous resistance. The wave plan is essentially a solution synthesized from 
the effects of incremental step pressure pulses. The pressure pulses are gen- 
erated because of incremental flow-rate changes that originate in a hydraulic 
system from a variety of sources, including the mechanical motion of the sys- 
tem structure. The pressure pulses propagate throughout the system at sonic 
velocity and are partially transmitted and reflected at each discontinuity. The 
velocity change caused by each pressure pulse is obtained from the Joukowski 
relation. Pressure and velocity time histories at any point in the system are 
obtained by a timewise summation of the contributions of the incremental pres- 
sure pulses passing that point. 

The analysis is presented in a form general enough to be applied to a 
variety of hydraulic systems. To illustrate the application of the method to a 
specific system, the response of a straight hydraulic line to a sinusoidal 
orifice-area variation of an upstream valve is computed. Both a constant- 
cross-section line and a tapered line are analyzed in the examples, and various 
nonlinear effects evaluated. Comparisons are carried out with experimental 
data obtained for the constant-diameter line and good agreement is shown to 
exist. 


INTRODUCTION 

Hydraulic systems employed in present day applications often require 
high performance with a small permissible range of variation from design 
flows and pressures. In addition, many hydraulic systems are structurally 
complicated. It is essential to be able to predict the dynamic response of these 
systems under transient and forced periodic conditions. 

*Asst. Prof, of Civil Engrg., Duke University, Durham, North 
Carolina. 
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The basic partial differential equations for unsteady flow are derived 
from continuity and momentum relations. Closed-form distributed parameter 
solutions for small sinusoidal flow perturbations in long lines having negligible 
fluid damping are available and are in good agreement with experimental 
•data*"*, In general, however, the nonlinear partial differential equation* can- 
not be solved without resorting to tedious numerical or graphic techniques. 

An analytical study was therefore undertaken at the NASA Lewis Re- 
search Center to develop a numerical distributed-parameter solution for un- 
steady flow in hydraulic systems that would be in a form readily lending itself 
to digital-computer computation. The analysis techniques which were devel- 
oped will be referred to herein as the wave-plan. The wave-plan method has 
some similarities to a method of characteristics solution due to the technique 
of tracing sonic disturbances throughout the systeih; however, the wave-plan 
method is more easily applied to complex unsteady flow systems in addition to 
lending itself better to physical interpretation. 

This paper presents the details of the wave-plan analysis along with the 
basic elements of the corresponding digital -computer program required to cal- 
culate unsteady flow in liquid systems. The analysis is presented in a form 
general enough to be applied to a variety of fluid systems. Examples are given 
to illustrate the application of the method to specific hydraulic systems. 


ANALYSIS 

General 


Basic equations . - The pressure and velocity of the liquid within a line as a function 
of position and time can be obtained from the equations of momentum and continuity. 

The momentum equation for a one-dimensional elastic fluid of constant mean density 


2S..l[22 + v2U /<v)| (i) 

3x g 3t 9x J 


where /(v) is the resistance due to fluid viscosity, which is some function of velocity. 
(Symbols are defined in appendix D. ) 

The continuity equation for an elastic liquid in an elastic line is^ 
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Equations (1) and (2) are nonlinear partial differential equations which, 
when solved along with the boundary conditions in a liquid-filled line, give re- 
lations for velocity and pressure in the line. 

Method of s olution. - The direct volution of the continuity and momentum 
Muatieni to Obtain vvTaattltfl And prvivurev is difficult even when the boundary 
conditions are relatively simple (open or closed ends). In addition, exact solu* 
tions have been found only for a system where the resistance term is linearized 
or neglected entirely^ - . 

Acceptable solutions of nonidealized unsteady flow problems including 
nonlinear effects require the use of numerical methods and a digital computer. 
The continuity and momentum equations have been solved by a numerical tech- : 
nique based on the method of characteristics and adapted for digital computer 
use^. Identical solutions may be obtained more readily by an alternate proc- : 
ess. This method, referred to as the wave-plan solution, is developed in this 
paper . 

A wave-plan solution is obtained as follows: At the point or points in a 
liquid system where a disturbance is introduced (such as an oscillating valve 
or a moving impedance change), an incremental change in liquid flow rate due 
to the disturbance over a very short interval of time is computed. The incre- 
mental pressure pulse accompanying the flow-rate change is then computed. 
This pressure pulse is propagated throughout the system at sonic speed. The 
pressure pulse is partially reflected and partially transmitted at all geometri- 
cal and physical discontinuities in the fluid network. Pressure and velocity 
time histories are computed for any point in the system by summing with time 
the contributions of incremental waves. 

Characteristic impedance . - The relation (characteristic impedance) be- 
tween pressure and velocity changes caused by a pulse traveling in the liquid- 
filled line is computed from momentum considerations. Figure 1 shows pres- 
sure and flow conditions, as a pressure wave propagates in a liquid-filled line, 
that exist at time t and at time t + At. The wave takes the time At to travel 
the distance Ax from A to B. During this time there is a pressure P + AP 
on the left end and a pressure P on the right end of the liquid contained be- 
tween points A and B. This unbalanced pressure causes the fluid to accelerate. 
Newton's second law gives 
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Figure L - Effect of pressure pulse on mean line conditions. 
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Canceling and rearranging yield 


AP 


p 

At 


The term Ax/ At la ths propagation apttd of the pressure wsvs, Tha wave apead la 
equal to the aonic velocity C in the ay stem if the mean velocity of the liquid in the line 
is neglected. Since the mean velocity of the liquid is usually much smaller than the 
acoustic velocity, this is usually permissible. Thus, 

AP - pC AV (3a) 

or in terms of head of liquid 


AH = — AV (3b) 

g 

The sonic speed C for a liquid flowing within a line is influenced by the elasticity of the 
line wall, and for a system that is axially unrestrained^**^ can be calculated 
from 



Equation (3), which is the well known Joukowski equation, states that the pressure 
change at any point in a line is equal to the product of the velocity change at the point 
times the characteristic acoustic impedance pC of the liquid in the line. This relation 

o o 

may also be derived by employing continuity and energy principle 


Generation of Pressure Waves 

Pressure waves may originate in liquid flow systems in a variety of ways. Two 
common sources are an orifice or valve with a varying open area and a moving point in a 
line where there is an impedance change. In order to apply the wave-plan method to 
varying area change or moving impedance change, the corresponding disturbance func- 
tion is approximated by a series of discrete changes over small equal time intervals. 
Figure 2 shows this approximation. 

Each of these changes is assumed to occur instantaneously at some time during the 
time interval. The pressure perturbations associated with each of these changes 
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are then computed as shown in the following examples. 

Varying area orifice. - A pressure wave is gen- 
erated at an orifice that experiences a change in open- 
ing area. The magnitude of the generated wave is a 
function of conditions before the area change and the 
magnitude of the am o hangs, 

The flow through an orifice is assumed to obey the 
square-law relation: 

V- B VaH 

The orifice coefficient B is a function of the open area, the line area, and the dis- 
charge coefficient and is easily determined for a particular orifice. For the case of an 
open area that is a small percentage of the line area (Aq « A), the orifice coefficient is 
given closely by 

B.c D ^yii 

A 

where the discharge coefficient Cq is a function of Reynolds Number (and hence, veloc- 
ity). For small velocity perturbations, however, the discharge coefficient varies only 
slightly, and hence the orifice coefficient may be considered to vary linearly with area 
changes. The examples in this report deal with this case. 

Figure 3 shows conditions at an orifice before and after a small orifice-area change. 
The flow is from left to right. The external pressure head may be varied in a pre- 
scribed manner. 

The momentum equation across the wave front gives 



Figure 2. - Step approximation of disturbing 
function. 


The flow out the orifice after the orifice area change is given by 
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Figure 3. - Effect of step change In orifice coeffi- 
cient al terminal orifice. 


V 11 - B 2 V »1 + AH 11 - «22 

c 

Solving the momentum and orifice equations simulta- 
neously gives the following quadratic in 

V 11 + bV ll + c “ 0 (8) 

where 
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The positive root of this equation gives the desired value for outflow velocity. Substi- 
tuting this value for V 11 back into the momentum equation gives the magnitude of the 
pressure wave generated by a sudden change in orifice coefficient from Bj to B 2 . 

Moving terminal orifice . - The analysis of the magnitude of a pressure wave gen- 
erated by the motion of a component in a fluid line is exemplified by considering the 
motion of a terminal orifice. The lateral velocity of the orifice is so represented by a 
series of step changes that the velocity of the orifice remains constant over short time 
intervals. Figure 4 shows conditions at a terminal orifice Just before and a short time 
after a step change in velocity. The orifice coefficient need not be constant. 

The momentum equation across the generated wave is 

AH n-f< v i- v n> 


The displacement shown in figure 4 takes place in the time increment At. The con- 
tinuity equation states that the net inflow across boundary AA equals the net flow out the 
orifice plus the storage that takes place during the time interval. Therefore 


V n A At = VEgA At + VOA At 


or 


v n = ve 2 + VO 
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Figure 4. - Effect of lateral motion of terminal orifice. 


where VE 2 is the lateral velocity of the ori- 
fice end during the time interval, and VO is 
the liquid velocity in the line (with respect to 
the orifice). 

The flow out the orifice is given by 
VO = B 2 yHj + AH n - H 22 

Solving the preceding equations simultaneously 
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for the resulting line velocity gives the following quadratic in VO: 


where 


VO 2 + bVO+c*0 



g 


and 


= -B^Ihj t SlVj - ve 2 ) - h 22 


Effect of Viscous Resistance 


After the generation of a pressure wave at a perturbing point in a liquid system, the 
wave propagates with sonic velocity throughout the system. The viscous resistance of 
the liquid medium influences the propagation of this wave. 

The effect of viscous resistance on pressure pulse propagation can be neglected 
without appreciable effect for short, large-diameter lines where such losses are small; 
however, for longer lines or small-diameter lines carrying viscous liquid, the resistive 
losses are not negligible and should be included in an unsteady flow analysis. 

The exact solution of the partial differential equations (1) and (2) is limited by the 
necessity of linearizing or neglecting the friction term entirely. For graphical analyses 
it has been necessary to consider the friction losses as lumped at one or more points in 
the pipeline if a manageable solution is to be obtained*®. 

To include the effects of viscous resistance in a wave-plan solution, the extent that 
fluid viscosity influences the propagation of a pressure pulse must be determined. 

There is presently no experimental or analytical information available which is in a 
form that can be used for the prediction of the effect of viscous resistance in unsteady 
line flow. Because of this limitation it is necessary to make a quasi-steady approxima- 
tion and to assume that viscous losses in a small line length dx are given at any instant 
by the Darcy equation as 


Ah 


_ f Ax V 
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( 7 ) 


The form of equation (7) is identical to that of an internal square-law orifice if the 
orifice coefficient is given by 
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This equation implies that viscous losses over a small length of line can be lumped 
at an internal square-law orifice with a properly chosen orifice coefficient. This repre- 
sentation is rsfsrrtd to as ths "or If to* analogy," and its use in graphical analysis has 
been suggested by Bergeron 10 . 

Line losses can be distributed at many discrete points along a line by inserting a 
large number of correctly chosen square-law (friction) orifices. Impinging waves are 
then reflected and transmitted at these friction orifices in a manner similar to that of 
reflection and transmission through a small region of flowing viscous liquid. The equa- 
tions governing wave reflection and transmission at a friction orifice are developed in 
the next section. 


Reflection and Transmission of Pressure Waves at System Discontinuities 


It is necessary at each discontinuity to compute the magnitude of transmission and 
reflection of each pressure pulse in terms of initial conditions at a discontinuity, the 
nature of the discontinuity, and the magnitude of the impinging waves. These deriva- 
tions are all made by employing the three fundamental fluid flow relations (continuity, 
momentum, and energy). 

Figure 5 shows two waves impinging and reflecting simultaneously at a discontinuity 
in a liquid system. The notation used in this figure is employed in subsequent specific 
examples. Liquid moving from left to right has a positive magnitude of velocity. 

The momentum equations across the wave fronts yield 



Figure 5. * Nomenclature for conditions before and alter wave 
action at discontinuity. 

(See fig. 5 for definition of V* and V M .) 


AH n = AHj 


AH i ' T <v ’ ‘ v i > 
iH n“T <v '- v ii ) 

AH 2 ' T <V 2 ■ V,,) 

a H 2 2 = T <V 22- V,,) 

Equations (8) and (9) combine to give 
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and equations (10) and (11) give 

ah 22 -ah 2 *3(v 22 -v 2 ) 
The new pressures are given by 


H 1X «Hj+ AHj + AHjj 

h 22 - h 2 + ah 2 + ah 22 


( 13 ) 


(14) 

(15) 


Terminal orifice . - A terminal orifice is considered to be bounded by a pressure 
reservoir. The pressure in the reservoir may be changing in a prescribed manner. In 
addition, the orifice coefficient may also be changing in a prescribed manner. Figure 6 
shows conditions before and after reflection from an orifice bounded at the inlet by a 
pressure reservoir. 
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Figure 6. - Pressure-pulse reflection at terminal 
orifice. 


If the discharge after wave action is in the posi- 
tive direction (from the reservoir to the line), the 
head-discharge relation for the orifice is given by 

V 22 - B 2 Vh u - (H, ♦ AHjs ♦ AH 22 ) 

Solving this equation simultaneously with the 
momentum equations (13) for the velocity gives the 
following quadratic in V 22 : 

V 22 + bV 22 + c = 0 (16) 


where 


and 


b = 


B 2 C 2 

g 


c = “ B 2 H 11 - H 2” 2AH 2 + 




Because the resulting direction of flow was taken as positive, the positive root of 
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equation (16) is the desired result. - 

The magnitude of the reflected pressure wave is obtained from equation (13). The 
new pressure at the orifice is given by equation (15). 

If the direction of flow is into the reservoir after wave action (negative direction), 
the velocity head relation is given by 

V„ - B, V<H a + AH, + AHjj) - H u 
Simultaneous solution with the momentum equations gives a quadratic in 

v 22 + by 22 + c * ® 

where 

g 


The negative root of equation (17) is the desired result. 

It is necessary to determine the resulting direction of flow so that the correct equa- 
tion (eq. (16) or (17)) can be applied. Inspection of equation (16) discloses that the term 
H 11 - H 2 “ 2 ah 2 + c 2 V 2^ g must 1)6 P° sitive to yield the necessary positive root. If 
this term is negative, the flow must be into the reservoir, equation (17) must be used, 
and the negative root must be computed. 

Friction orifice . - Figure 7 shows conditions before and after pulse-wave action at 
a friction orifice. Because of identical conduits on both sides of the orifice, the head- 
velocity relation for the orifice after wave action for flow from left to right is 

v„ - V 22 . B f yHj + AHj + AH n - Hjj - AH 2 - AH 22 
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Figure 7. - Pressure-pulse reflection it internal orifice. 


Noting that = C 2 = C and solving this equation 
simultaneously with the momentum equations (12) 
and (13) give 

V 11 + bV ll + c * 0 (18) 

where 
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2B 2 C 
b = — L_ 

g 

2 / 2CV i 

c - -BpfHj + 2 AHj - Hj - 2 Attj + * 

The positive root of this quadratic equation gives the resulting velocity throu gh the 
friction orifice. The term Hj + 2 AH^ - H 2 - 2 AH 2 + 2CVj/g must be positive to yield 
a positive root. If this term is negative, the resulting flow is from right to left, and the 
head-velocity relation for the orifice is 

V U • V 22 * Bp V H 2 * 4H 2 + iH 22 ' H 1 ‘ iH l ‘ 4H 11 

This equation combines with the momentum equations to give the following: 

V ll + bV ll + c “° (19) 



where the coefficients b and c are of the same magnitude and of opposite sign from 
the corresponding coefficients in equation (18). 

The negative root of this equation gives the resulting velocity in the negative direc- 
tion. The magnitudes of the resulting pressure waves and the pressures after wave 
action are given by equations (12) to (15). 

Diameter discontinuity . - The relation between net head and discharge for a diam- 
eter change in a conduit is derived by utilizing energy relations. Figure 8 shows the 
notation for this derivation. 

If A j > Ag, the diameter discontinuity is an abrupt constriction, and the energy 
equation for flow through the constriction from left to right after wave action is 
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Figure 1 - Pressure-pulse reflection at diameter discontinuity. 
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The area ratio is denoted as 



Then the continuity equation is 


v 2 * 


( 20 ) 


11 



V 22“ RV 11 


( 21 ) 


If these relations are substituted into the energy equation and the terms are com- 
bined, the following net head-velocity relation is obtained: 




( 22 ) 


If A ^ < Ag, the diameter discontinuity is an abrupt expansion. The energy equation 
for flow through the expansion from left to right after wave action is 


w2 tr2 

V H „ v 22 .. 

— — * + H. . ** — + Ha a + 

2g 11 2g 22 
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This equation gives the following velocity-head relation for an expansion: 


H U- H 22= V 




( 23 ) 


The case of a lossless diameter discontinuity is also important because step diam- 
eter changes can be employed to simulate a tapered line. In the continuously tapered 
line the losses are small and sometimes can be neglected. The energy equation without 
the loss term is 

V 2 V 2 

V ll +H - 22 „ 

IF » IT 22 


The resulting velocity-head relation is 


H 11 - H 22 ' V ll( 5 17 1 ) 


( 24 ) 


The final net head after wave action can be written as 


H n - H 22 = Hj + AHj + AH U - H 2 - AH 2 - AH 22 


or 
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H n - H J2 - Hj + 2 AH, + Y (V 1 ' V U> - H 2 - J 4H 2 - ^ (V 22 - V 2 ) (25) 


The relation is combined with the continuity equation and the velocity-head relation 
to obtain second-order polynomials for the line velocity after wave action* The equation 
i> _ 


aVu + bVjj + c ■ 0 


(26) 


where 
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The coefficient a depends on the type of discontinuity and whether energy losses 
are included. For a contraction uith losses included the coefficient is 


a = 


3R* - 2R - 1 
4g 


(27) 


For an expansion with losses included the coefficient is 


„ _ R- R 
a = 

g 


(28) 


For either type of discontinuity with no losses included the coefficient is 


a = 


R 2 - 1 
2g 


(29) 


For each case the resulting velocity is the positive root of equation (26). There- 
fore, the term 


, c 2 v 2 


Hj + 2 AHj + ^ + 


- Ho - 2 AH, 


g g 
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must be positive. If this term is negative, the resulting flow is from right to left, and 
the equations are reformulated by putting the loss term on the correct side of the energy 
equation. 

Junctions of three or more legs . - The previously presented mathematical analyses 
pertaining to the reflection of a pressure pulse at a diameter discontinuity have taken 
into aooount energy loesee introduced by the discontinuity. With Junctions of three or 
more legs, the inclusion of energy loss terms becomes Increasingly difficult because of 
the irreversibility of frictional effects. Thus, the treatment employed herein will be 
similar to the one presented in standard references on fluid dynamics and will not be 
derived**. 

The relations for computing the percentages of magnitude of a pressure wave trans- 
mitted and reflected at a junction are based on the following: 

(1) The Joukowski equation applied across each wave 

(2) Continuity of flow at the junction 

(3) Continuity of pressure at the junction 

The magnitude of the wave that is transmitted to all the other legs of a junction of 
n legs due to a wave of magnitude AH impinging in line i is given by T(i)AH, where 
the transmission coefficient T(i) is given by 

T(i) = - (30) 



j=l 


The magnitude of the wave reflected in leg i is given by R(i)AH when the reflec- 
tion coefficient R(i) is given by 


R(i) =T(i) - 1 


Analysis of the Dynamics of a Liquid System by Wave Plan 

The wave-plan analysis supposes a system composed of a discrete number of dis- 
continuities connected by lossless line segments, which serve only to transmit pressure 
pulses. The discontinuities include geometric ones, such as terminal orifices and diam- 
eter discontinuities, and artificial ones, such as friction orifices. A typical simulated 
fluid system is shown in figure *9. 
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Thtfi 1« a vuliWi-un orlftst at A (valve) and a flxad orlflee at B. At C there is 

an abrupt change in line diameter. Friction orifices are inserted between A and C and 
between C and B to simulate viscous resistance. 

The equations developed thus far have applied only to conditions at a particular point 
in a system at a particular instant of time. Through the use of these equations it is 
possible to compute velocity, pressures, ami magnitudes of pressure pulses leaving a 
particular discontinuity in terms of the magnitudes of impinging waves and conditions 
^rior to wave action. 

To analyze a system, however, the equations must be both time and position de- 
pendent. This dependency is indicated by the form of the basic partial differential equa- 
tions (eqs. (1) and (2)). It is apparent that the pressure wave impinging at one discon- 
tinuity in a liquid system is exactly that wave which left an adjacent discontinuity at some 
time in the past (because of the lossless line connector simulation). Thus, the waves 
are related by position and time. In general, variables have the form /(x,t). Specifi- 
cally, H(x,t), V(x,t), and AH(x, t) denote the pressure head, velocity, and pressure 
pulse as a function of position and time, with x and t being the position and time sub- 
scripts, respectively. 

The wave -plan analysis provides for information at specified points in the system 
(discontinuities) and at discrete time intervals. New information is available only at 
discrete time intervals because all of the disturbing functions are approximated by a 
series of small step changes occurring at specified times. The system response to 
these disturbing functions will also be in the form of step changes. Because of the two- 
parameter dependence of the dependent variables, it is advantageous to denote position 
and time subscripts. 

Position subscripts . - In order to apply the analysis equations that have been devel- 
oped for a discrete point to a liquid system, it is necessary to introduce a method of 
denoting position. This is done by numbering adjacent discontinuities in a consecutive 
manner as indicated in figure 9. (This arrangement is modified at junctions of three or 
more legs). The position subscript is denoted as L. 

In addition, because pressure and velocity states may differ across discontinuities, 
it is necessary to denote a left and a right side of a discontinuity. This is done by adding 
an L for left and an R for right to the nonsubscripted part of the variable. 

Thus HL(L, t) denotes the pressure head to the left of discontinuity L at time t. 
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Time subscripts . - The time subscript J is defined as follows: 

t - J At 

where t is the time and At is the working time increment. Thus HR(L, J) denotes the 

pressure head to the right of discontinuity L at time t ■ J At. 

Selection of working time increment . - The working time increment At is the time 
interval between successive computations. 

Its selection is determined by two factors. First, the time increment must be 
small enough to approximate accurately all disturbing functions by a series of step 
changes. Second, it is necessary that all reflections of pressure pulses in the system 
take place at an integer number of time intervals. The wave travel times between all 
adjacent discontinuities will then be an integer number of working time increments. If 
this were not the case, waves would be impinging on a discontinuity in a completely 
arbitrary fashion and would make the solution unmanageable. 

The selection of the time increment is simplified because the majority of discon- 
tinuities in a system will be friction orifices, and these can be placed where desired. 

For example, the selection of the working time interval for the system shown in figure 9 
is made as follows: 

(1) The wave travel times between adjacent geometric discontinuities are determined: 

(a) tac - wave travel time between A and C 

(b) tcb - wave travel time between C and B 

Within the desired limits, the largest time interval of which these travel times are 
integer multiples is determined. This time interval represents the largest working time 
increment possible. 

(2) It is then determined if the integer number of working time increments necessary 
for travel between A and C and between C and B is large enough to insert the desired 
number of friction orifices that must, of course, be separated in time by at least one 
working time increment. If not, the increment can be divided by any integer to obtain a 
smaller working time increment. 

(3) It must be determined if the working time interval necessary to assure an integer 
number of time intervals between all discontinuities is small enough to approximate 
accurately the disturbing function (variable-area orifice at A) by a series of step changes. 
If not, the working time increment is further divided by an integer necessary to get a 
suitable approximation. The integer number of time increments between discontinuities 
is increased by this factor. 

Subscripted notation for analysis equations . - The analysis equations are easily ex- 
tended for application to a liquid system by the use of the subscripted notation. Fig- 
ure 10 shows conditions at a discontinuity before and after wave action when the sub- 
scripted notation is employed. 
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Figure XL - Subscripted nomenclature for conditions before and alter wave action at 
discontinuity. 


The time that the waves reach discontinuity L is t = J At. Conditions before and 
after that time are constant, and step changes occur at t « J At. 

A comparison of the nonsubscripted notation in figure 5 to the subscripted notation 
in figure 10 gives the identities necessary to make the analysis equations general expres 
sion 8 for a liquid system. These identities are as follows: 

VL(L, J - 1) = Vj = velocity to left of discontinuity L 
HL(L, J - 1) = = pressure head to left of discontinuity L 

VR(L, J - l) s v 2 s velocity to right of discontinuity L 
HR(L, J - 1) * H 2 « pressure head to right of discontinuity L 
AHR(L - 1, J - IOC) = AHj = pressure pulse coming from adjacent discontinuity to 
left. KX is the number of working time intervals it takes a sonic disturbance 
to travel between the two discontinuities. 

AHL(L + 1, J - KY) = AH 2 = pressure pulse coming from adjacent discontinuity to 
right. KY is the number of working time intervals it takes a sonic disturbance 
to travel between these two discontinuities. 

The conditions after wave action are as follows: 

VL(L, J) = Vjj = velocity to left of discontinuity L 

HL(L, J) = = pressure head to left of discontinuity L 

VR(L, J) = V 22 - velocity to right of discontinuity L 

HR(L, J) = H 22 = pressure head to right of discontinuity L 

AHL(L, J) = AHj j - pressure pulse leaving discontinuity L and moving to left 

AHR(L, J) = AH 22 = pressure pulse leaving discontinuity L and moving to right 

Substitution of these identities into the analysis equations (eqs. (8) to (29)) yields a 
perfectly general set of algebraic equations applicable to any liquid system. 

Method of computation . - The solution is carried out by using the subscripted equa- 
tions as follows: 

(1) A working time increment is selected. 

(2) Initial conditions are computed from given eteady- state values that exist in the 
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system prior to the initiation of a disturbance. These conditions include velocity and 
pressure to the left and right of each discontinuity and are denoted at t = 0 by VL(L,0), 
VR(L,0), HL(L,0), and HR(L,0). 

(3) Computations are then carried out to determine conditions at each discontinuity 
at the end of the first working time Interval by using the conditions from Btep (2) as ini- 
tial conditions. It should be noted that the analysis equations will predict no change at an 
undisturbed discontinuity that is not subjected to impinging waves. Thus, conditions will 
only change at the end of the first time interval at a disturbance source (point A, fig. 9). 
At the disturbance source conditions change, and waves are generated that will subse- 
quently alter conditions at adjacent discontinuities. 

(4) Each discontinuity is analyzed at the end of each interval. The computations are 
carried out as long as desired to give pressures and velocities at each discontinuity at 
the end of each working time interval. 


Digital Computer Programing 

The wave -plan analysis consists of the sequential solution of many equations. The 
dynamic analysis of even a rather simple liquid system for a reasonable number of time 
intervals involves a large number of computations. A digital computer must be employed 
if the solution is to be obtained in a reasonable length of time. p 

To make the digital analysis as general as possible, computer segment programs 
(subroutines) have been written for the most common types of discontinuities. These 
discontinuity solutions can then be combined (incorporating junction equations, if neces- 
sary) to obtain digital computer models for various liquid flow systems. 

Each discontinuity subroutine computes conditions at the discontinuity for some 
point in time as a function of the local velocity and pressure-head conditions a time in- 
terval earlier, the magnitude of impinging waves (from adjacent discontinuities), and the 
step changes in disturbing functions during the time interval. 

The following subroutines have been written to be used in the digital model repre- 
sentation: 

(1) Terminal orifice subroutine 

(2) Internal (friction) orifice subroutine 

(3) Diameter discontinuity subroutine 

Each subroutine uses the general relations given by equations (12) to (15). Also 
each subroutine initially assigns nonsubscripted identities to their subscripted counter- 
parts and finally assigns the nonsubscripted results to the subscripted terms (as indi- 
cated in the section entitled Subscripted notation for analysis equations). A description 
and computer printouts of the digital computer subroutines are included in appendix A. 
All digital computer programming has been done in Fortran IV. 
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The subroutines are combined to form digital computer models for the dynamic 
analysis of liquid flow systems. Two examples are included in the next section. 

RESULTS AND DISCUSSION 

The equations and programs developed in the ANALYSIS section form a basis for the 
analysis of unsteady flow conditions in a liquid flow network. The equations and pro- 
grams have been kept general so that they may be applied to the analysis of a variety of 
liquid systems. In general, the wave-plan digital-computer analysis has distinct advan- 
tages over other methods such as analog- computer models that substantially linearize or 
lump parameters and small perturbation techniques that linearize around mean-line con- 
ditions. The following features of the wave-plan analysis should be noted: 

*1 (1) Viscous friction effects are easily included. 

(2) Perturbing functions of any form may be used. 

(3) Nonlinear relations are easily included. (Linearization of parameters is of little 
or no advantage. ) 

(4) Both transient and steady-state responses to disturbing functions are given. 

(5) Complex networks can be analyzed. 

Two examples will be used to illustrate the application of the wave-plan method to 
the construction of digital models of particular liquid flow systems. The first example 
will illustrate the method of computing the transient and steady-state response of a long, 
straight hydraulic line to a periodic sinusoidal flow disturbance. In the second example 
the response of a long tapered line will be computed. 


Example 1: Long Hydraulic Line With Oscillating Inlet Orifice 

This example was chosen because the computer results can be compared 
directly with experimental data from the Lewis Research Center line dynamics 
r ig4-6. Schematics of the liquid flow system are given in figure 11. Fig- 
ures 11(a) and (b) show the analytical model and the experimental rig, respec- 
tively. The system consists of a 68-foot-long, 7 /8-inch-inside-diameter line 
with pressure reservoirs at both ends. The test fluid was JP-4 fuel or another 
hydrocarbon, depending on the liquid viscosity desired. Line pressures 
ranged from 200 to 400 pounds per square inch. A variable-area throttle valve 
that perturbs the system is located at the upstream end of the line (point A), 
and a fixed orifice is located at the downstream end (point B). The details of 
constructing the digital model are given in appendix B along with the program 
constants used in the calculations. 
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Inlet t- Orifice with Outlet 



Figure II. - Hydraulic line with oscillating Inlet orifice. 

Effect of amplitude of input perturbations . - Most line-dynamics analyses are based 
on small-perturbation techniques that permit the use of linearized terminal impedances 
as an approximation for nonlinear pressure-flow relations at the end of the line. There- 
fore, in taking experimental line-dynamics data^~ ^ the practice is to limit the 
amplitude of the oscillating throttle disturbance generator. The minimum usable ampli- 
tude, however, for a given line condition and terminal impedance is limited in practice 
by the ratio of the sine wave signal to the nonharmonic noise which decreases with de- 
creasing amplitude. The permissible maximum amplitude is limited in practice by the 
appearance of harmonics in the sinusoidal pressure signal as the amplitude is increased. 
The system must be operated between these two limits. 

The relation between the disturbance amplitude and the degree of nonlinearity of the 
pressure and velocity perturbations was determined analytically by varying the amplitude 
of the sinusoidal input perturbations for the analytical line model. The mean orifice 
coefficients were 0. 8 at the inlet and 1. 0 at the outlet. The steady line velocity was 
14 feet per second. For the digital computer program (appendix B) the amplitude of the 
input orifice coefficient perturbation BA was varied at 40 cps as shown in table I. 
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The calculated steady-state responses for the 
pressure and fluid velocity at the inlet and outlet of 
the line (points A and B, fig. 11(a)) are shown in 
figure 12 (p. 22). The velocities and pressures are 
dimensionless, having been divided by steady line 
value*. Figure 18 show* that, m the amplitude of 
the disturbing function is decreased, the system 
responses are more nearly described by a 
sinusoidal-type periodic function; however, for the 
particular line condition and terminal impedance 
used in the example, nonlinearities are evident 
even for relatively small perturbations. For ex- 
ample, the inlet velocity perturbation is very clearly nonsinusoidal for an 18. 75-percent 
disturbing amplitude even though the perturbation amplitude is only 2 percent in the posi- 
tive direction and 4 percent in the negative direction. 

Transient response . - Dynamic systems analyses are usually based on steady-state 
responses to steady sinusoidal inputs. This analysis technique does not, however, pre- 
dict the transient response to the initiation of a disturbance or a change in the disturb- 
ance characteristics in a liquid flow system. The wave-plan solution gives both trans- 
ient and steady-state solutions. This point is clarified by examining both the transient 
and final steady-state responses to the initiation of a sinusoidal disturbing function in a 
liquid system under steady mean-flow conditions. A value for BA of 0. 15 (case 4, table I) 
was chosen for the amplitude of the disturbance function. The line response starting with 
the initiation of the disturbing function is given in figure 13 (pp. 23 and 24). For the case 
studied several cycles were required for the perturbation velocities and pressures to ap- 
proach the steady-state values. In some cases the amplitudes of the transient responses 
are considerably different from the steady responses. For example, during the early 
part of the transient period the magnitude of the inlet velocity perturbation was nearly 
double the steady -state value. 

Comparison of analytical and experimental results . - The experimental hydraulic 
line (fig. 11(b)) was operated with the upstream throttle valve sinusoidally perturbed at 
70 cps with area perturbation amplitudes of 12. 3, 33. 8, 55. 3, and 76. 8 percent of the 
open area. The mean value of the coefficient for the upstream orifice was 0. 65, and the 
corresponding value for the downstream orifice was 0. 42. The downstream orifice was 
slightly open with respect to characteristic. Figure 14 (p. 25) shows the resulting ex- 
perimental pressure-time oscilloscope traces for the quarter-length point of the line 
(measured from the upstream throttle valve). 

The experimental conditions were duplicated by the digital model, and figure 15 
(p. 25) shows the theoretical results for these four different input amplitudes. The 


TABLE I. - AMPLITUDES OF ORIFICE 
COEFFICIENT PERTURBATIONS 


Case 

Input orifice 
coefficient 
perturbation 

amplitude, 

BA 

Percent of steady- 
state orifice 
coefficient 

1 

0.48 

60 

2 

.32 

40 

3 

.20 

25 

4 

.15 

18. 75 


\ 
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Figure 12. - Effect of amplitude of Inlet orifice area perturbation on steady-state response of long line. 

agreement between the wave shapes determined analytically and experimentally is very 
good for all cases. In addition the analytical and experimental magnitudes agree to 
within the accuracy of the measurements (approx. 3 percent). 

This example illustrates that the wave-plan analysis for unsteady liquid flow is 
capable of accurately including nonlinear effects and giving results in good agreement 
with experimentally observed values. 
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Example 2: Response of Long Tapered Line to Sinusoidal input Disturbance 


The analysis of the response of a liquid system composed of lines of nonuniform 
diameter (tapered) to a periodic disturbing function is difficult if classical methods are 
employed; however, tapered lines can easily be included in a wave-plan analysis by 
approximating the tapered line by a line composed of a discrete number of diameter 
changes (fig. 16, p. 26). 

In order to evaluate the quantitative effect of tapered lines, the system shown in 
figure 17 (p. 26) was studied in detail. This system consists of two pressure reservoirs 
connected by a tapered line. At the end of the tapered lines, relatively large resistive 
losses were introduced in the form of square-law orifices. The throttle valve at point A 
has a sinusoidally varying area, and the orifice at point B is fixed. 
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Figure 13 . - Concluded. Transient response of long line to Initiation of sinusoidal Inlet valve area pertur- 
bation. 
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Figure 14. - Experimental pressure-time oscilloscope traces for quarter-length point of line. 




Time 

Or) Orifice area perturbation amplitude, 33. 8 percent of open 

area. 




(d) Orifice area perturbation amplitude, 76. 8 percent of open 
area. 


Figure 15. - Analytical pressure-time traces for quarter-length point of line. 
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Flfur* 16 . • Stap-dtfmtfer-^tenft igrotMtten to ttptrtd lint. 


This system was studied for various 
degrees of taper. The configurations 
were chosen to give the diameter at the 
inlet as 

DA-D • 6 ft (»U) 


and the diameter at the outlet as 


DB = D + 6 ft (31b) 

The diameter of the tapered conduit 
is assumed to vary linearly between the 
ends, and the average diameter is always 
D feet. 

If the resistances of the terminal orifices are large compared to the line loss, the 
steady discharge of the system is practically independent of 6 and would depend mainly 
on the reservoir pressures and the mean line diameter. Even for lines with much 
smaller end-resistive losses the system discharge would be only slightly dependent on 
the degree of taper, because line losses consisting of friction and expansion or contrac- 
tion losses would be mainly dependent on line length and average diameter. Mathemati- 
cally this can be written as 



q2 * HEN- HEX (32) 

K + K, + K 
en l ex 

where K gn , Kp and K gx are the loss coefficients of the entrance orifice, the line, and 
the exit orifice, respectively. 

When this system is used with large values of K en and K gx as compared to Kp 
it is possible to compare directly the dynamic response of systems that are essentially 
statically equivalent. Pressure and flow perturbations can be quantitively compared as 
a function of degree of taper. The straight conduit (6 = 0) is a special case that is in- 
cluded in the analyses. 

The tapered conduit was chosen so that the wave velocity would be constant; thus the 
ratio of diameter to wall thickness remained constant along the length of the line. 

For the digital model the tapered line was approximated by N - 1 short straight 
sections of conduit, as shown in figure 17. In this manner the diameter of each approxi- 
mating section I is computed by the equation 
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The details of the digital model along with the computer flow diagram 
and program constants are given in appendix C. One basic system was studied 
in detail, a system having a high-resistance orifice (three-fourths to seven- 
eights closed geometrically but slightly open with respect to characteristic im- 
pedance) at the inlet with a 10-percent sinusoidal variation in the area of the 
orifice opening. The outlet orifice was chosen of even higher resistance than 
the inlet and represents an end that is closed with respect to characteristic 
impedance. Nine degrees of taper were studied, the average diameter D 
being 1 foot for each case. The values of 6 that were used were -0.3, -0.2, 
-0.1, 0, 0.1, 0.2, 0.3, 0.4, and 0.5 foot. 

The quarter-wave resonance frequency for the nontapered line is 15 cps. Frequen- 
cies of 5, 7. 5, 10, 12. 5, 15, 17. 5, 20, and 25 cps were studied for each case. In all 
cases the pressure and flow perturbations were slightly nonlinear even though they were 
small. In all cases they were within 2 percent of mean-line values. Figure 18 shows 
typical pressure and flow perturbations for one cycle. 
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Id! Pressure-perturbation gain tor contracting taper. 
Figure 19. - Frequency dependence of taper ed-llne dynamic response. 


In spite of the nonlinearity of the perturbations, they were steady and repeatable 
after several cycles (during which the transient died out) and resembled sinusoids. It 
was considered appropriate to make comparisons based on positive amplitudes. This 
procedure is followed in the ensuing discussion. 

The analytical results are displayed graphically in figure 19. Figure 19(a) shows 
the variation of the amplitude of the outlet pressure perturbation with frequency for the 
expanding taper (line that expands from inlet to outlet). Figure 19(b) shows the pressure 
perturbation gain (ratio of outlet pressure perturbation amplitude to inlet pressure per- 
turbation amplitude) for the expanding taper. Figures 19(c) and (d) show the same re- 
sults for contracting tapers. The most striking result apparent in figure 19 is the shift 
in the resonant frequency for tapered lines. This shift is shown in figure 20. 

Another result is that the maximum pressure perturbations at the outlet have a 
strong dependence on 6 and vary considerably over the entire frequency range. Thus a 
certain amount of control can be exerted over the response of this system by the judi- 
cious choice of tapers. For example, in the case of an expanding taper of 6 = 0. 3 the 
amplitude of the outlet pressure perturbation is significantly reduced when compared to 
a straight-line response for the frequencies shown except in the range from 6 to 12 cps. 
A line with a contracting taper of 6 = -0. 3 differs little from a nontapered line (6 = 0) 
up to 15 cps. For frequencies greater than 15 cps the outlet pressure amplitude pertur- 
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Figure 20. - Variation ot resonant frequency wttti taper factor. 


bation is considerably greater. 

This example indicates that the 
dynamic response of a liquid system 
can be altered significantly with a cer- 
tain amount of control without signifi- 
cantly altering the static response of 
the system; however, the primary 
intent of the example is to illustrate 
the use of the wave-plan analysis. 
Although the example was applied to 
a system with a linear taper and a 
sinusoidal disturbing function, a sys- 
tem of any arbitrary taper with any 
disturbing function could have been 
analyzed with no additional difficulty. 


CONCLUSIONS 

The analytical method presented in this paper provides a means of ob- 
taining distributed parameter solutions to a variety of unsteady flow problems 
for liquids flowing in conduits. An advantage of the method is that the com- 
plete solution is obtained. For example, both the transient and the steady- 
state response to a suddenly imposed periodic flow disturbance is obtained. 
Furthermore, the disturbing function can be of arbitrary form and need not be 
periodic. Nonlinear effects are easily included. In addition, the wave-plan 
method is advantageous in making certain types of dynamic response calcula- 
tions. For example, the response of liquid-filled lines having types of axial 
cross-sectional area distributions for which there would be little hope of ob- 
taining closed-form analytical solutions can be easily handled. 

The wave-plan analysis can readily solve problems in which there is 
interaction between the structural motion of the conduits and the perturbations 
in the fluid flowing within the conduits. Interactions of this type often occur 
in hydraulic systems; however, the motion is usually neglected in the analyses 
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APPENDIX A 
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COMPUTER SUBROUTINES 

Digital '■computer routines have been written to solve for conditions at various line 
discontinuities after wave action in terms of the conditions at the discontinuity prior to 
wave action, the magnitudes of impinging waves, and the physical characteristics of the 
discontinuity. These routines have been formulated in such a manner that they may be 
easily incorporated into a computer program to synthesize different liquid flow systems. 

In all cases the computer routines have been written in nonsubscripted notation, 
which corresponds exactly to the notation presented in the text of this paper. In the 
calling vector that calls the subroutine into the main program, the subscripted counter- 
parts of the nonsubscripted subroutine variables are identified. 

This mechanism of identifying subscripted variables in a calling vector 
makes it possible to use the analyses presented in the text for specific orienta- 
tions and flow conditions to solve for all orientations and flow conditions. In 
this manner the analyses presented in the text of this paper have been general- 
ized. Proper identification of the variables in the calling vector is tantamount 
to orienting the discontinuity to correspond to the case analyzed in the text. A 
table incorporated into each subroutine gives the proper subscripted-to- 
nonsubscripted dependence for each case. 

Each of the analyses presented in the text resulted in a polynomial that 
is quadratic in a velocity term. For several reasons these equations were 
solved by an iterative manner by employing Newtonian extrapolation. The roots 
of the polynomials were determined by looping the following equations: 


aVll 2 + bVll + c 

error = 

2aVll + b 
Vll = Vll - error 

The looping is continued until the error is of sufficiently small magnitude. An ac- 
ceptable solution can be obtained within a very few cycles because the wave-plan solution 
deals with small changes, and the value of the velocity before wave action is a good 
approximation of the final value. In addition, the coefficients b and c are generally 
large and always much larger than a. These conditions lead to a rapidly diminishing 
error term. 

A direct solution of the quadratic equations by the quadratic formula results in a 
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solution that is a small difference between two large numbers. In soiqe cases sufficient 
accuracy cannot be obtained even when double-precision computing techniques are em- 
ployed. The computer programs which use notation identical to that appearing 
in the text of this paper didactically illustrate the logic of the subroutines and 
are presented in lieu of flow diagrams. ' A brief description of the program and 
the Fortran IV program for each computer routine follows . 

Terminal Orifice Subroutine 


Equations (16) and (17) are the basic equations employed in this analysis. They are 
derived in the text for an orifice, on the left end of a line. The table at the beginning of 
the subroutine gives the identities that must be substituted into the calling vector depend- 
ing on whether the orifice is on the right or left end of the line. Before the subroutine 
is entered, the following terms must be identified with a numerical value: 

B orifice coefficient 

HEN reservoir head for left end 
HEX reservoir head for right end 

J time counter 

K number of time increments to nearest discontinuity 

L position counter 

The computer program for a terminal orifice is as follows: 


SUBROUTINE T0R(H2 »DH2 *V2 ,H11 *H22 »DH22 .V22 ) 
SUBROUTINE TERMINAL ORFICE 
COMMON/BOX/A, B,C» G*BF .ERROR 


c 


LEFT END 

RIGHT END 

c 

H2 

HR(L.J-l) 

HL ( L » J-l ) 

c 

V2 

VIL.J-1 ) 

-VIL.J-1 ) 

c 

DH2 

DHL ( L+l *J-<) 

DHRIL-l.J- 

c 

V22 

V ( L * J ) 

-V ( L » J ) 

c 

H 22 

HR(L*J) 

HL ( L » J ) 

c 

DH22 

DHR(L.J) 

DHL (L*J) 

c 

HI 1 

HEN 

HEX 


BB = 

C*B*B/G 



CC=-B*B* ( HI 1-H2-2 • *DH2+C*V2/G I 


IF(CC) 1 0 * 10 * A 
A BB = -BB 

cc = -cc 

10 V22 = V2 

11 ER = (V22*V22+BB*V22+CC)/(2.*V22+BB) 
V2 2 = V22-ER 

IF ( ABS ( ER ) -ERROR ) 12*12.11 

12 DH22=DH2+C/G* < V22-V2 ) 

H22=H2+DH2+DH22 

RETURN 

END 
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Friction Orifice Subroutine 


The basic equations for this routine are equations (18) and (19). Before this routine 
is entered, the following terms must have a numerical value: 

BF friction orifice coefficient 

J time counter 

K number of time increments to nearest discontinuity 

L position counter 

The computer routine is the following: 

SUBROUTINE FOR ( Hi >DHl »DH2 *H2 »V1 »Hl 1 »DH22 »DHl 1 *H22 *V11 ) 

C SUBROUTINE FRICTION ORFICE 


c 

HI 

HULtJ-l ) 

c 

DH1 

DHR (L-l » J-K ) 

c 

DH2 

DHL(L41,J-K) 

c 

H2 

HR ( L » J-l ) 

c 

VI 

V ( L * J— 1 ) 

c 

HI 1 

HL t L * J ) 

c 

DH 1 1 

DHL I L ♦ J ) 

c 

DH22 

DHR (L ♦ J ) 

c 

H22 

HR(L»J) 

c 

Vll 

VIL.J) 


C0MM0N/B0X3/ A * B* C*G*BF»ERRDR 
BB * 2.*BF*BF*C/G 

CC*-( Hl+2 »*DHl+2 »*C/G*Vl-( H2+2 »*DH2 ) )*BF*BF 
IFtCC) 3»3*2 

2 BB * -BB 
CC * -CC 

3 Vll * VI 

4 ER * (Vll*Vll+BB*Vll+CC)/tBB+2.*Vll) 

Vll « Vll-ER 

IF (ABSIER)-ERROR ) 5»5»4 

5 DH1 1 * DH1+C/GM Vl-Vll ) 

HI 1 * H1+0H1+DH1 1 

DH22 * DH2+C /G* ( V 1 1— VI ) 

H22 * H2+DH2+DH22 

RETURN 

END 

Diameter Discontinuity Subroutine 


Equations (20) to (29) are employed in the diameter discontinuity subroutine. These 
equations are derived for flow in the positive direction (to the right) for both an expan- 
sion and a contraction. Before the subroutine is entered, the direction of flow after 
wave action is determined, and this direction in turn determines which subscripted var- 
iables are inserted into the calling vector. 

The subroutine also includes the case where no viscous loss is considered across 
the discontinuity. To utilize this case, it is necessary to assign the term LOSS a value 
of zero. Additional terms that must be defined are the following: 


32 


nnnnnnnnnn nAastt^ non 


J time counter 

K number of time increments to nearest discontinuity 
L position counter 
The computer routine is as follows: 


tUHQUTINI D !IC ( Hi tVl * PHI »Hi»VliSH!»Ci *C1*A1»AI »Hll *DHi i »Hll »DHlI » 
1V11.V22) 

C DISCONTINUITY SUBROUTINE 


COMMQN/BOX/A »B »C * G*BF *ERR0R 
C0MM0N/B0X2/ LOSS 



FLOW TO RIGHT 

FLOW TO LEFT 


CCC GREATER THAN 0 

CCC LESS THAN 

HI 

HL ( A* J-l ) 

HR ( L ♦ J-l ) 

VI 

VLIL.J-1 ) 

-VR(L*J-1> 

DH1 

DHR ( L-l * J-K > 

DHL ( L+l ♦ J-K ) 

H2 

HR(L»J-1 ) 

HL<L»J-1> 

V2 

VR (L . J-l ) 

VL ( L » J-l > 

DH2 

DHL(L+1 *J-K> 

DHR(L-1»J-K) 

HI 1 

HLCL.J) 

HR ( L * J ) 

Vll 

VL ( L * J ) 

-VR(L»J) 

DH1 1 

DHL ( L » J) 

DHR(L.J) 

H22 

HR ( L » J ) 

HL(L.J) 

V22 

VR t L » J ) 

-VL(L.J) 

DH22 

DHR( L*J) 

DHL(L.J) 

Cl 

CL ( L ) 

CR(L) 

C2 

CR ( L ) 

CL(L) 

A1 

AL t L ) 

AR (L) 

A2 

AR(L) 

AL(L) 

R = A1 

/A2 


CC* 

“(Hl+2.*DHl+Cl*Vl/G-H2-2. 

*DH2+C2*V2/G> 

BB * 

C1/G+C2*R/G 


IF (LOSS) 4*3.4 


AA = 

(R*R-1. )/<2.*G) 


GO TO 10 


I F ( R 

-1.) 1*1*2 


AA * 

(R*R-R) /G 


GO TO 10 


AA * 

(3.*R*R-2.*R-1. >/(4.*G> 


Vll* 

VI 


ER = 

( AA*V11*V11+BB*V11+CC)/(2.*AA*V11+BB) 

Vll* 

Vll-ER 


IF ( ABS( ER ) -ERROR ) 6*6*5 


DH1 1 

* DH1+C1 /G* ( Vl-Vl 1 ) 


V22* 

R*V1 1 


DH22 

*DH2+C2/G* ( V22-V2 ) 


Hi 1 = 

Hl+DHl+DHl 1 


H22* 

H2+DH2+DH22 



RETURN 

END 


3 
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APPENDIX B 


DIGITAL MODEL OF HYDRAULIC LINE WITH OSCILLATING INPUT ORIFICE 

A schematic drawing of the system chosen for this example is shown in figure U 

(p. 20). 

The system consists of a variable-area orifice at A that perturbs the system and a 
fixed orifice at B. The conduit is bounded by pressure reservoirs. (The reservoir 
pressures can vary. ) 

To keep the solution general, N - 2 friction orifices are inserted, which result in 
N discontinuities. A time interval is so chosen that each discontinuity is one time in- 
terval from adjacent discontinuities. Thus 

Atf.W-U 

C \N - 1/ 

• This time increment is a multiple of the working time increment At, which is 
chosen small enough so that all disturbing functions can be accurately described by step 
functions. 

Initial values for velocity through each discontinuity and the pressure to the right 
and left of each discontinuity must be inserted into the computer as data or the computa- 
tion of these quantities from the initial steady condition must be integrated into the com- 
puter program. 

The disturbing function can either be read into the computer as a function of time or 
it can be computed. For this case a sinusoidal orifice coefficient for the inlet (at A) is 
taken to be 


B = BO + BA SIN(2ir FJAt) 

The orifice coefficient for the outlet is considered to be constant and is denoted 
by BN. 

The following data are needed as input to the computer: 

C wave velocity in conduit, ft/sec 

L length of line, ft 

BO mean orifice coefficient when subjected to sinusoidal perturbations 
BOl steady-state inlet orifice coefficient 

BA amplitude of inlet orifice coefficient perturbation 

F frequency of pressure perturbation 
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BN outlet orifice coefficient 

VO steady line velocity, ft/sec 

HEN pressure head of inlet reservoir, ft 

HEX pressure head of outlet reservoir, ft 

2 

g acceleration due to gravity, ft/sec 

N number of discontinuities 

NN number of working time increments for wave to travel length of conduit 

K number of working time increments between discontinuities 

M total number of time increments for which computations will be made 

N2 number of friction orifices for which pressure and velocity data will be 

printed 

In setting up the computer program the notation J * x, y, z is used. The notation 
means that computations are carried out for j starting at j = x. The computations are 
repeated for values of j that are increased by y until j > z. 

The following output parameters are printed out by the computer: 

V(1,J), HR(1,J) 

V(N, J), HL(N, J) 

V(N2, J), HR(N2, J), HL(N2,J) 

For a specific example, the fluid system constants appearing in the computer pro- 
gram were chosen nominally to match experimental conditions in the Lewis line- 
dynamics facility shown schematically in figure 11(b) (p. 20). Two sets of line condi- 
tions were used. The first (case I) 

table n. - constants for example was chosen in conjunction with an 

analytical study of the effect of the 
size of the input amplitude on the 
linearity of the perturbation imposed 
on the mean flow and pressure in the 
line. The second (case n) was chosen 
to match conditions frequently obtained 
in the line. 

The constants listed in table II 
were employed. A simplified 
computer flow diagram for this 
example is presented in figure 21. 
A complete listing of the Fortran 
IV program for this example is 
available^. 



Case I 

Case n 

C, ft/sec 

3800 

3600 

L, ft 

68 

68 

BOl * BO 

0.8 

0. 65 

BA 

0. 48,0.32,0. 20,0, 15 

0. 08,0. 22, 0. 36,0. SO 

BN 

1,0 

0. 42 

VO, ft/sec 

14 

6.5 

HEN, ft 

520 

360 

HEX 

0 

0 

g, ft/sec 2 

32.2 

32.2 

N 

17 

17 

K 

1 

2 

M 

170 

400 

NN 

16 

32 

F, cps 

40 

70 

N2 

5 

5 

ERROR 

0.00001 

0. 00001 
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APPENDIX C 


DIGITAL MODEL FOR TAPERED HYDRAULIC LINE 

The system chosen for the tapered hydraulic line if shown In figure 17 (p. 26). This 
system consists of a variable-area orifice at A and a fixed orifice at B connected by a 
tapered line with a linear diameter variation between A and B. The line is bounded by 
pressure reservoirs. The tapered line was approximated by a line composed of a dis- 
crete number of constant -diameter sections of decreasing (or increasing) diameter. 

To keep the solution general, N - 2 diameter discontinuities are inserted, which 
■result in N discontinuities and N - 1 sections of straight conduit. 

( The conduit is so divided that the wave travel times between adjacent discontinuities 
are multiples of the working time interval, and the diameter of the approximating 
straight section of conduit is equal to the average diameter of the section that it is 
approximating. 

It is assumed that the wave velocity remains constant over the entire conduit, which 
is equivalent to assuming that the ratio of the conduit diameter to wall thickness remains 
constant (a reasonable assumption from a structural viewpoint). The conduit is thus 
approximated by N - 1 sections of equal length. The diameter of the approximating 
straight section is given by 

/fa B - Da\/ 2I - l\ 

D(D=DA + (- 1 — j 

Initial values for velocity and pressure to the right and left of each discontinuity 
must be inserted into the computer, or the equations for the computation of these quan- 
tities from initial steady conditions must be integrated into the computer program. 

The disturbing function can either be read into the computer as a function of time 
or it can be computed. For this model a sinusoidal orifice coefficient for the inlet (at A) 
is taken to be 


B = BO + BA SIN(2vFJAt) 

The orifice coefficient for the outlet is constant and denoted by BN. 
The following data are needed as input to the computer: 

DA diameter at A, ft 
DB diameter at B, ft 
CA wave velocity in conduit, ft/eec 
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L length of line, ft 

BOl steady-state input orifice coefficient 

BO mean orifice coefficient when subjected to sinusoidal perturbations 
BA amplitude of orifice coefficient perturbation 

F frequency of pressure perturbation 

BN output orifice coefficient 

HEN pressure' of input reservolri ft 

QO steady-state flow rate, ft 

g acceleration due to gravity, ft/sec 

N number of discontinuities 

NN number of working time increments for wave to travel length of conduit 
K number of working time increments between discontinuities 

M total number of time increments for which computations will be made 
N2 number of diameter discontinuity for which pressure and velocity data will 
be printed 

The following output parameters are printed out by the computer: 

VRU,J), HR(1, J) 

VL^J), HL(N, J) 

VL(N2,J), HL(N2, J) 

HR(N2, J) 

For a specific example the effect of line taper on the response of a particular sys- 
tem of the type shown in figure 17 (p. 26) was computed. Calculations were made for 
nine different degrees of line taper. The specific data used in this analysis were as 


follows: 


TABLE m. - END DIAMETERS AND 
CONDITIONS FOR TAPER STUDIES 


DA 

DB 

BO 

BA 

BN 

1 

1 

1 

0.1 

0.2 

.9 

1.1 

1.235 

.1235 

.165 

.8 

1.2 

1.562 

.1562 

.139 

.7 

1.3 

2.041 

.2041 

.1185 

.6 

1.4 

2. 78 

.278 

.102 

.5 

1.5 

4.0 

.4 

.089 

1.1 

.9 

.825 

.0825 

.247 

1.2 

.8 

.695 

.0695 

.313 

1.3 

.7 

.593 

.0593 

.408 


CA = 3000 ft/sec 
L = 50 ft 
HEN = 660 ft 
QO =6.4 cu ft/sec 
N = 11 
NN = 10 
M = 300 

Table m gives the end diameters and conditions 
for the various tapers studied. 

These cases were studied for frequencies of 
5, 7. 5, 10, 12. 5, 15, 17. 5, 20, 22. 5, and 25 cps. 

Figure 22 gives a simplified computer 
flow diagram for this example. A listing 
of the Fortran IV program for this example 
is available 1 ^. 
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APPENDIX D 




SYMBOLS 


A 

lint arta, aq ft 

A(i) 

line area of 1, sq ft 

A(j) 

line area of j, sq ft 

A 0 

orifice area, sq ft 

A 1 

line area to left of diameter dis- 
continuity, sq ft 

a 2 

line area to right of diameter 


discontinuity, sq ft — 

a,b,c 

polynomial coefficients 

B 

orifice coefficient, (ft/sec 2 ) 1 / 2 

b f 

friction orifice coefficient, 
(ft/sec 2 ) 1/2 

B 1 

orifice coefficient before wave 
action, (ft/ sec 2 ) 1 / 2 

B 2 

orifice coefficient after wave 
action, (ft/sec 2 ) 1 / 2 

c 

sonic velocity in line, ft/sec 

C D 

orifice discharge coefficient 

C(i) 

wave velocity in line i, ft/sec 

C(j) 

wave velocity in line j, ft/sec 

C 1 

wave velocity to left of discon- 
tinuity, ft/sec 

C 2 

wave velocity to right of discon- 
tinuity, ft/sec 

D 

line diameter, ft 

D 

mean diameter of tapered line, 
ft 


S e , clastic modulus of oonvoyor, 

lb/sq ft 

&£ elastic modulus of fluid, lb/sq ft 

f Darcy friction factor 

/ function 

g acceleration due to gravity, 

ft/sec 2 

H pressure head, ft 

HEN pressure head of inlet reservoir, 

ft 

HEX pressure head of outlet reser- 
voir, ft 

Hj pressure head to left of discon- 
tinuity before wave action, ft 

H 2 pressure head to right of discon- 
tinuity before wave action, ft 

**11 pressure head to left of discon- 

tinuity after wave action, ft 

H 2 2 pressure head to right of discon- 

tinuity after wave action, ft 

AH pressure head change across 

wave or orifice, ft 

AHj wave impinging from left of dis- 

continuity before wave action, 
ft 

AH 2 wave impinging from right of 

discontinuity before wave 
action, ft 
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AH 11 

wave leaving discontinuity on 

V 

velocity in line, ft/sec 

left side after wave action, ft 

, VE X 

average verity of moving end 

AH 22 

wave leaving discontinuity on 


before wave action, ft/sec 


right side after wave action, 

VE a 

average velocity of moving end 


ft 


after wave action, ft/aec 

Ah L 

pressure head loss ever small 
line length, ft 

VO 

velocity In line adjacent to mov- 
ing orifice, relative to orifice, 

J 

time subscript 


ft/sec 

K 

number of working time incre- 
ments between discontinuities 

V 1 

velocity in line to left of discon- 
tinuity before wave action, 

L 

line length, ft 


ft/sec 

P 

pressure in line, lb/sq ft 

V 2 

velocity in line to right of discon- 
tinuity before wave action, 

AP 

pressure wave, lb/sq ft 


ft/sec 

Q 

volume flow rate, cu ft/sec 

V 11 

velocity in line to left of discon- 

R 

area ratio at diameter discon- 
tinuity 

tinuity after wave action, 
ft/sec 

R(i) 

reflection coefficient in line i 

V 22 

velocity in line to right of discon- 

T(i) 

transmission coefficient in line i 

tinuity after wave action, 
ft/ sec 

t 

At 

t‘ 

time, sec 

incremental time interval, sec 
short time before wave action, 

X 

6 

position, ft 

line taper factor (eqs. (31)), ft 


sec 

P 

mass density of fluid, slugs/ft^ 

t + 

short time after wave action, sec 

r 

wall thickness of conduit, ft 
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“ / Th rough A 

m/ ~ 


Conditions at 
left end 
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B - BO + BA $in(2nFJ At) 


h0 


K Terminal / Through \ f Friction 

orifice B i \ or ' fice 

subroutine^ / \l - 2, 1, N - 1/ \subroutine 
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orifice 
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Figure 21. - Computer flow diagram for hydraulic line with oscillating input orifice. 


start 




Figure 22. - Computer flow diagram for tapered hydraulic line with oscillating input orifice. 








